Langevin Dynamics simulations of a 2-dimensional colloidal crystal under confinement 

and shear 
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Langevin Dynamics simulations are used to study the effect of shear on a two-dimensional colloidal 
crystal confined by structured parallel walls. When walls are sheared very slowly, only two or three 
crystalline layers next to the walls move along with them, while the inner layers of the crystal are 
only slightly tilted. At higher shear velocities, this inner part of the crystal breaks into several 
pieces with different orientations. The velocity profile across the slit is reminiscent of shear-banding 
in flowing soft materials, where liquid and solid regions coexist; the difference, however, is that in 
the latter case the solid regions are glassy while here they are crystalline. At even higher shear 
velocities, the effect of the shearing becomes smaller again. Also the effective temperature near the 
walls (deduced from the velocity distributions of the particles) decreases again when the wall velocity 
gets very large. When the walls are placed closer together, thereby introducing a misfit, a structure 
containing a soliton staircase arises in simulations without shear. Introducing shear increases the 
disorder in these systems until no solitons are visible any more. Instead, similar structures like in 
the case without misfit result. At high shear rates, configurations where the incommensurability of 
the crystalline structure is compensated by the creation of holes become relevant. 



I. MOTIVATION 

Colloidal particles are a very useful model system for 
the study of fundamental properties of condensed matter, 
for example for the study of disorder and statistical fluc- 
tuations [l|-13|. These particles are convenient, because 
their interaction is tunable in many ways. For instance, 
by adding variable concentrations of non-adsorbing poly- 
mers to a suspension of colloidal particles, their hard-core 
repulsion is amended by a depletion attraction whose 
strength is controlled by the polymer concentration. An- 
other example are charged colloids, where adding salt can 
modify the strength and range of the effective Yukawa 
potential between the particles. Finally, if one uses col- 
loidal particles that contain a superparamagnetic core, 
a magnetic field can be used to control the strength of 
the r~'^-dipole-dipole-potential acting between them [l|- 
Q. Their size also makes them useful for experiments 
as it typically lies in a range which can be observed di- 
rectly by confocal miscroscopy. Additionally, the size and 
shape of these particles can be varied. It is also feasible 
to trap them in layers \Vm. [9l-[l3j. e.g. at air- water in- 
terfaces, in order to get a two-dimensional system. Such 
systems can be confined in a mechanical way or by laser 
fields d [13, Q 111- The effect of shear on the homoge- 
neous and heterogeneous nucleation of colloidal suspen- 
sions [3]3], on phase separation [20] and on colloidal 
glasses [21[ has also been adressed. 

One of the important questions in this context is how 
the structure formation in colloidal systems can be in- 
fiuenced by the various tunable parameters (interaction 
potential, size, shape, mixtures of different kinds of par- 
ticles) and by imposing external conditions like confine- 



ment, shear, compression and temperature [22]. A thor- 
ough understanding of these processes is necessary for the 
design of methods to create nanostructures utilizing the 
self-assembly of particles. The aim here is to be able to 
predict which pattern will form under which conditions. 

Several studies using Monte Carlo computer simula- 
tions have already been carried out investigating the 
structure formation in colloidal crystals. Some of them 
have focussed on the effect of different kinds of walls 
|23l [24| , on the impact of external fields 25] and on con- 
finement incommensurate with the crystalline structure 
leading to the formation of a soliton staircase (26l - t28| in 
a one-component colloidal system crystallizing in a tri- 
angular lattice structure. Additional studies have inves- 
tigated the melting [2^, the deformation [30], the cor- 
relation between channel geometry and pattern forma- 
tion [Slj, [3^ and frustration [S^ of such two-dimensional 
systems of purely repulsive particles under confinement. 
Others have considered binary mixtures which crystallize 
into a square lattice structure [34| - |40| and fiuid systems 
[41| . In this paper we will concentrate on the influence of 
shear on a one-component colloidal crystal with a hexag- 
onal lattice structure with and without misfit. The aim 
of our work is to clarify the range over which moving 
boundaries induce fiow of particles in such a soft crystal, 
and to study the extent to which such nonequilibrium 
phenomena can still be described by local equilibrium 
concepts. 

In the following section, we will describe the model sys- 
tem and the simulation parameters. In section Hill we will 
report on the results of our simulations with sheared walls 
for the case where the effective periodic wall potential is 
commensurate with the lattice period of the crystal. We 



will present snapshots of the created configurations, ve- 
locity profiles and angular distributions. Then we will 
move to the results for the case where the lattice period 
is no longer commensurate with the wall potential and 
present the same quantities in order to compare them. 
As a last step, we will look at the mean square displace- 
ment of the particles in both cases. In section ITVl we will 
briefly summarize the main conclusions of this work. 



II. DETAILS OF THE SIMULATION 

In order to investigate the influence of shear in com- 
bination with confinement in a two-dimensional colloidal 
crystal, we set up a model system consisting of iV = 
73440 particles, arranged in a slit geometry which is 
elongated along the x-axis. In the x-direction, periodic 
boundary conditions are applied, while the system is con- 
fined in the y-direction by two walls consisting of two 
rows of frozen particles each (see fig. [T] These "frozen" 
particles interact with the "real" particles via the same 
potential with which the particles interact with each 
other. We use 30 rows parallel to the walls (in a com- 
mensurate arrangement with the walls) or only 29 rows 
(which then are incommensurate due to higher particle 
densities in the resulting smaller "volume" ) . 
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FIG. 1. Sketch of the system geometry with the directions 
of the axes and the sheared wall particles. Black dots denote 
wall particles, grey dots freely movable particles. 

Shear is implemented by moving these structured walls 
in the x-direction with a certain, constant velocity. The 
two walls on either side of the stripe are moved with the 
same velocity, but in opposite directions. 

The particle interaction is defined via a purely repul- 
sive pair potential like the potential which would be felt 
by magnetic particles when a magnetic field is turned on. 
However, the resulting potential is inconvenient for com- 
puter simulations due to its slow r~^ decay with distance, 
and since there are not yet any experiments for sheared 
colloid crystals to compare with, we prefer to use a com- 
putationally more efficient r~^^-potential with a cutoff 
(in order to make it strictly short-ranged) , a shift (in or- 
der to ensure that it ends with a value of zero), and a 



smoothing factor (to make it differentiable) . A further 
motivation to use this model potential is that it is essen- 
tially the same potential used in the studies of confine- 
ment effects on colloidal crystals without shear |23l - [28| 
and of two-dimensional melting (42] . This potential is 
defined as M,M: 
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where h — O.Olcr and a = 1 were chosen. Thus the 
particle diameter a is our unit of length, and e — I sets 
the scale for the temperature. 

The particles are arranged in a hexagonal lattice struc- 
ture. In the case without misfit, the walls are placed 
such a distance apart, that a perfect crystalline lattice 
structure with ny — 30 rows of n^ — 2160 particles fits 
between them. The dimensions of the system L^ and 
D follow from the lattice constant a via L^ = n^a and 
D = nj,av3/2. In the case with a misfit, the walls are 
placed closer together, thus shrinking the width of the 
stripe to D — {uy — A)a\/3/2. This parameter A is re- 
ferred to as "misfit" in the following. 

At larger misfits, the structure of the crystalline strip 
is rearranged to 29 rows as has been investigated in pre- 
vious studies by Chui et al [H-ii]. At a misfit of A = 2.2 
we therefore started a part of our simulations with a com- 
pressed 30-row-structure, which quickly rearranged itself 
into a 29-row-structure upon shearing, and a part of our 
simulations were started with a 29-row-structure. This 
structure with 29 rows contained as many particles as the 
structure with 30 rows. Therefore, those particles that 
had been in the 30th row had to be distributed amongst 
the remaining 29 rows. As the number of particles in the 
walls did not change, previous studies showed, that the 
rows directly adjacent to the walls remain free of extra 
particles while all of the inner rows incorporate the same 
number of extra particles thus creating a lattice with a 
smaller lattice constant ax in x-direction. But due to the 
misfit, the lattice constant in y-direction is smaller, too, 
which stabilized this new crystalline structure. Never- 
theless, the incommensurability of the crystalline struc- 
tures of the inner rows with the rows directly adjacent 
to the walls leads to the formation of a soliton staircase 
[26j . Fig. [5] explains how solitons are created by the in- 
commensurability of two rows with a different number of 
particles. 

The simulations were carried out using the program 
package HOOMD-blue 5^, |4J1' which is designed to run 
on graphic cards and provides the option to use Langevin 
Dynamics by only defining a parameter 7 for the drag of 
the particles, a timestep At and a temperature T. It 
uses the well-known Velocity Verlet integration method 
and adds a force F — —7-? -|- Fraud to the force exerted 
on each particle by the interaction with its neighbouring 
particles [4^, \^. Here, v is the particle's velocity and 
Fraud is a random force with a magnitude chosen via the 
fluctuation-dissipation theorem to be consistent with the 



FIG. 2. Explanation of solitons: As there are more particles 
in the the inner rows, the distance between these particles is 
on average smaller than the distance between the particles in 
the outer rows forming the rigid walls. These rigid rows act 
like a given external periodic potential on the mobile particles. 
Therefore the particles in the inner rows are no longer placed 
at the minima of the potential created by the particles in the 
outer rows, thus creating a soliton. 



specified drag 7 and temperature T. In all of the simu- 
lations presented here, we used 7 = 0.5. Note that our 
particles have mass m = 1, and time t is measured in 
the standard Molecular Dynamics time unit t = i^/— ^. 
Because of the huge speed-up that occurs when MD sim- 
ulations are run on graphic cards we initially used them 
as well, but experienced problems with the accumula- 
tion of systematic errors. However, as the graphic card 
code uses single instead of double precision, we had to 
carry out the simulations on regular CPUs as the higher 
accuracy turned out to be necessary for our crystalline 
systems. We chose a timestep of At = 0.002 and several 
temperatures as indicated in the respective paragraphs. 
All of these temperatures were below the melting tran- 
sition which would occur for this system at the chosen 
density oi p = 1.05 (for the case without misfit) at about 
T — 1.35 [2aj. Nevertheless, the Langevin thermostat 
only takes into account the velocities of the particles 
in the bulk (and not the velocity of the wall particles 
which are moved with constant velocity in order to in- 
troduce shear into the system). This leads to a higher 
effective temperature of the particles due to their inter- 
action with the moving wall particles, because the wall 
particles are constantly introducing friction and thereby 
increasing the temperature (see fig. [3]). The Langevin dy- 
namics are aiming at restalling the initial temperature, 
but as the wall particles move again with every new time 
step, more heat is added to the system, and an equilib- 
rium between the cooling due to the Langevin thermostat 
and the heating due to the sheared walls is reached. This 
equilibrium temperature is higher than the temperature 
which was set in the Langevin thermostat. But this tem- 
perature increase due to the sheared walls does not only 
depend on the velocity of the walls, but also on the inter- 
actions of the particles with the walls. So at rather low 
shear velocities between v = 0.1 and v = 3.5, the effective 
temperature of the system increases with increasing shear 
velocity, but for a very high shear velocity oi v — 50.0, 
the effect of the shearing on the temperature is almost 
negligibly small again, because such a fast shearing has 
almost no effect on the particles next to it, which feel 
only an averaged-out wall potential in this case. 

Since near the walls the effective temperature ex- 
tracted from the x and y-components of the velocities 
of the particles disagree, it is clear that the "fluid lay- 
ers" near the moving walls are completely out of equilib- 
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FIG. 3. Effective temperature, calculated from the velocity 
distribution in x-direction (a) and y-direction (b) by fitting 

it via f{vx) = ao ■ e'"""""""' '' "ff' (and likewise for Vy, with 
Vy =0). Shear velocities of the walls as indicated. Starting 
configuration with 30 rows without misfit (A = 0), at T = 
0.3. The inset in part (a) compares the effective temperatures 
Teff.x (solid lines) and Teff,y (dotted lines) for 11 — 5.0, v = 
1.0, V = 0.5 and v = 0.1 (from top to bottom). 



rium, since the concept of a "local temperature" which 
is anisotropic does not make much sense. However, for 
small velocities the "temperature" profiles T^ff^xiv) and 
Ti,ff^y{y) agree within statistical errors, and the con- 
cept of a "local temperature" in the slit interior can be 
used. Then the observed profile of this local tempera- 
ture Te//(y) can be interpreted as follows: due to the 
large friction between the fiuid-like layers and the mov- 
ing walls, these layers act like a source where heat is 
constantly pumped into the system. This heat is trans- 
ported via conduction into the interior of the slit (due 
to the action of the thermostat, heat is removed every- 
where in the slit, and a steady-state temperature profile 
is established). Motivated by corresponding solutions of 
the heat conduction equation, we can fit these local tem- 
perature profiles by T{y) = AQ[exp{— z / X) + exp{—{D — 
z)/X)] + Al, where Ao, A\ and A are phenomenological 
parameters that depend on v. Unfortunately, due to the 
large fiuctuations the fitting of three parameters is some- 
what uncertain and hence we do not dwell on the velocity 
dependence of these parameters. It would be interesting 



to study how this dependence changes when the friction 
constant is varied, but a study of this problem would 
reqiure very large computational resources and hence is 
beyond the scope of the present paper. 



III. RESULTS 

Without a misfit (i.e. A = 0), where the crystal con- 
sists of 30 rows, slow shearing of the walls has a rela- 
tively large effect on the velocities of the particles. Fig.|3] 
shows velocity profiles obtained by sorting the particles 
into rows depending on their position in the y-direction 
and then calculating the average particle velocity for each 
row of particles. As the maximum average particle ve- 
locity is given by the shear velocity of the walls, we nor- 
malized the velocity profiles by the shear velocity. But 
we also show the non-normalized velocities in an inset of 
the graphs (fig. |4]) as a higher absolute velocity will lead 
to more disorder in the crystal, while the normalized ve- 
locities indicate how strong the effect of the shearing is. 
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FIG. 4. Velocity profiles. Shear velocities of the walls as 
indicated. Starting configuration with 30 rows without misfit 
(A = 0), for two temperatures T = 0.3 (a) and T = 0.5 (b). 

At a temperature of T = 0.3 and for the case of A = 
this results in velocity profiles (fig. |4^) where the two 
smallest shear velocities {v = 0.1 and v — 0.3) lead to 
the (relatively) largest particle velocity in the first 4 — 5 
rows closest to the walls. This means that slow shearing 



of the walls has the strongest effect on pulling the par- 
ticles along with the walls. Of course, if one compares 
the absolute, non-normalized velocities inside the crystal, 
these velocities are larger at a shear velocity of w = 0.3 
and V — 0.5, as particles next to the walls obtain a higher 
velocity due to the faster walls, even if a smaller percent- 
age of the walls' velocity is transferred to the particles 
next to them. The higher absolute velocity also leads 
to more disorder inside the crystal, which then leads to a 
stronger influence of the shearing even on particles in the 
inner rows of the crystal. Thus, there is clear evidence 
for shear band formation, as will be discussed below. 

This disorder is also reflected in the density profiles 
(fig. EI and in the angular distributions (fig. [7^). The 
density profiles show clear peaks indicating a crystalline 
structure for higher shear rates, but at a shear rate of 
V = 0.3, the peaks corresponding to the inner rows of 
the crystal are almost non-existent any more. This can 
be explained by the changes in the crystalline structure 
discussed below. 
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FIG. 5. Density profiles. Simulations were carried out at a 
temperature of T = 0.3 for a system with 30 rows without 
a misfit (A = 0). Shear velocities of the walls as indicated. 
The strongly oscillating profile represents « = 0.8 and was 
omitted in the inset, where a magnified plot of the central 
part of the slit is shown. The normalization of the density 
profile p{y) is given in terms of the total number of particles, 

Obviously the profiles of the velocity never resemble 
the (approximately) linear velocity profiles representing 
simple Couette flow, that one encounters when one shears 
simple [131 or complex [H, Hi] fluids. This is due to the 
fact that we are shearing a crystalline solid here. At 
higher temperatures and small shear velocities almost 
linear velocity profiles can be obtained for our system 
as well (see fig.[9l). 

The shearing of the walls leads (at certain shear rates) 
to a tilting of the crystalline layers and even to breaks 
in the crystalline structure, where parts of the crystal 
are oriented in a different direction from the rest. This 
effect is illustrated in the snapshots shown in fig. HI In 
order to analyze it quantitatively, we computed for each 
pair of two neighbouring particles the angle between the 
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FIG. 6. Snapshots of a section of a system with sheared walls. 
The starting configuration contained 30 rows and the simula- 
tion was carried out at a temperature T = 0.3 and without 
a misfit (A = 0). The particles in the two first rows on top 
and the two last rows at the bottom (that are at rigidly fixed 
positions relative to one another and contribute the moving 
walls) are highlighted by black dots while the mobile parti- 
cles are shown as grey dots. Shear velocities of the walls were 
chosen as u = O.f (a), v = 0.3 (b) and v = 0.8 (c) 



line connecting these two particles and the x-axis. In 
a perfect hexagonal lattice without misfits and without 
defects, this would lead to peaks at 0°, 60°, 120° and 
180° (we did not compute any angles larger than 180°). 

At a shear velocity of w = 0.1, the peaks are centered 
around slightly larger angles due to a tilting of the crys- 
talline layers which is induced by the shearing (fig.|S^ and 
fig. [Zl). At shear velocities oi v — 0.3, v = 0.5 and (less 
pronounced) oi v = 0.8 the angular distribution has two 
peaks as the crystal breaks into several pieces and gener- 
ally only the layers close to the walls are aligned with the 
walls (creating a peak in the angular distribution close 
to 60°, 120° etc), while the crystalline parts inside of 
the crystal prefer a different orientation, thus creating 
the second peak in the angular distribution (fig. [6)d and 
fig.jT]). These two peaks are very distinct at ti = 0.3, while 
at u = 0.5 there is more disorder inside of the crystal due 
to the higher shear velocity. Starting from v = 0.8, this 
second peak disappears again as the walls are moving too 
fast to have much influence on the crystalline structure 
any more. This is due to the fact that the fast shearing 
of the walls leads to an averaged out wall potential which 
the particles next to the wall feel, which is not able to 
drag particles along with the walls any more. 

It is interesting to note that the above results may 
be understood as follows. Firstly, at low velocities, the 
velocity profile shown in fig. |3] may be interpreted as be- 
ing composed of three prominent shear bands |50| , a fea- 
ture of a typical non-Newtonian, shear-thinning material. 
The two fast moving bands lie close to the boundaries and 
fiow in opposite directions sandwiching a central immo- 
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FIG. 7. Distributions of angles between the line connecting 
each pair of two neighbouring particles and the x-axis. Shear 
velocities of the walls as indicated. The graphs for different 
shear velocities are shifted vertically by 4000 in order to give a 
better overview. Starting configuration with 30 rows without 
misfit (A = 0), temperature T = 0.3 (a) and T = 0.5 (b). 



bile region. In some respects our flowing solid therefore 
behaves like a moving paste, or soft-glassy matter [5l| 
in a fashion which is complementary to that seen in the 
experiments of Coussot et al. [52l | where a central mo- 
bile region is flanked on either side by immobile shear 
bands adjacent to the boundaries. As the shear rate is 
increased, the width of the mobile shear bands initially 
increases as expected. In order to verify that we indeed 
do have shear banding in our system, we plot the average 
stress as a function of shear rate in flg. [8] The resulting 
flow curve shows non-monotonic behavior typical of shear 
thinning systems. However, surprisingly, instead of ob- 
taining a linear Newtonian region at high shear rates, 
the stress stabilizes to a constant signifying the presence 
of a shear band of vanishingly small viscosity. Referring 
back to the calculated velocity profiles, (see fig. S]) we 
observe that the width of the shear bands decreases with 
increasing velocity in this regime, so that at large veloci- 
ties (shear rate), we obtain thin, highly mobile boundary 
layers near the walls while most of the solid is immobile. 
For such velocities, the solid re-crystallizes. 

Further, in order to characterize the behavior of the 
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FIG. 8. Flow curve for T = 0.3 and starting configurations 



with 29 rows at a misfit of A = 2.2. a^ 
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note the respective components of the stress tensor, while v 
indicates the velocity with which the walls are being moved. 
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FIG. 9. Velocity profiles for T = 1.5. Shear velocities of 
the walls as indicated. Starting configuration with 30 rows 
without misfit (A — 0). 



solid at high wall velocites, we refer to the plot of the 
anisotropic, local, effective temperature (obtained by fit- 
ting a Gaussian to the probability distribution of the 
X and y components of the velocity) already shown in 
fig. [31 We observe that both these effective tempera- 
tures show non-nionotonic behavior, initially increasing 
and then decreasing as the velocity increases leading to 
re-crystallization. Regaining crystalline order with high 
drive velocities is also seen in crystals driven over ran- 
dom pinning potentials [53|] . This implies that our system 
never shows a Newtonian regime at any fiow velocity if 
the temperature is small enough to obtain crystallization. 
At higher temperatures, when the solid melts in equilib- 
rium, we do obtain (nearly) Newtonian flow (fig. [9|). 

At higher temperatures, but still below the melting 
temperature, it is energetically easier for the crystal to 
split up into several domains, which can then be tilted 
with respect to the walls and can be oriented in different 
directions. Therefore the velocity profiles at temperature 
T — O.b have a more pronounced shape even at higher 
shear velocities (see fig. HJa) . The same effect is visible in 



the angular distributions (fig. [T)^) : The "second peak" in 
the angular distributions only disappears at considerably 
higher shear rates than in the case of T = 0.3. Apart from 
this, all shear rates at T = 0.5 lead to higher average 
particle velocities than in the case of T = 0.3 as the 
higher temperature makes it easier for particles to move 
along with the walls. 

A different picture arises if one introduces a misfit into 
these simulations. As described above, a misfit means 
that the structured walls are placed closer together. In 
simulations without shearin g, th is would lead to the cre- 
ation of a soliton staircase [2g (compare fig. [2]). It is 
also already known that there is a significant hysteresis 
in the transition between structures with 30 rows and 
structures with 29 rows with the same overall number of 
particles. These structures with 29 rows are the ener- 
getically preferred configuration for misfits greater than 
about A = 1.7, but if a simulation is started with a 
structure of 30 rows, it will not spontaneously change 
its structure to 29 rows (in simulations without shear) 
until about a misfit of A = 2.0. But even at larger mis- 
fits, simulations without shear have to be equilibrated 
for a certain time until this transition takes place. This 
transition happens much faster when sheared walls are 
introduced into the system. For example, for a system 
starting with 30 rows at a misfit of A = 2.2 and at a 
temperature of T = 0.3, the system still consisted of 30 
rows after 1 mio. time steps in the case without sheared 
walls, but had already 29 rows after just 20000 time steps 
in the case of a simulation where the walls were sheared 
with a velocity of u = 0.5. As will be discussed later on, 
at this shear velocity the structure breaks and changes 
again if the simulation run is continued. 

In the following, we will concentrate on simulations at 
a misfit of A = 2.2, which leads to a structure with 29 
rows independently of whether the starting configuration 
contained 30 or 29 rows. It is interesting to investigate 
the structures which develop here as they do not contain 
solitons any more due to the shearing of the walls. In- 
stead, it is now possible for the system to re-arrange its 
particles continually in a way that no particles have to 
sit on the maxima of the potentials, as there is enough 
disorder in the system, especially close to the walls. 

There are several mechanisms of dealing with the dif- 
ferent number of particles in the different rows which lead 
to the formation of solitons in the case without shear. At 
low shear rates and if the starting configuration already 
contains 29 rows, the row next to the walls (which con- 
tains the same number of particles as the wall rows, but 
less than all of the inner rows), moves along with the 
walls, while the inner rows remain almost without move- 
ment (fig.Hni curves for v = 0.1 and v = 0.3). This would 
in principle lead to solitons which are moved along with 
the walls as the rows directly adjacent to the walls still 
do not have the same crystalline structure like the inner 
rows. But as these rows are moving, the particles they 
contain are not always at the ideal lattice positions and 
therefore the solitons are strongly "smeared out" . 
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FIG. 10. Velocity profiles. Shear velocities of the walls as 
indicated. Temperature T — 0.3. Starting configuration with 
29 rows at a misfit of A = 2.2. 



If one starts with a structure of 30 rows at the same 
shear rate, the structure changes into 29 rows. Apart 
from equilibration effects, we obtained the same results 
as in the case where we started with a well-defined 29- 
row-structure. 
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FIG. 11. Angular distributions. Shear velocities of the walls 
as indicated. The graphs for different shear velocities are 
shifted vertically by 4000 in order to give a better overview. 
Starting configuration with 29 rows at a misfit of A = 2.2. 
Temperature T — 0.3. 



At slightly higher shear velocities (v — 0.5 and, less 
pronounced, also v — 0.8), the same effect as in the 
case without misfit is observed: The crystalline layers 
are tilted and the inner part of the crystal breaks off 
from the rest and changes its orientation, thus creating 
a rather strong flow in the layers close to the wall and 
less flow inside of those parts of the crystal, which have 
turned around. In this case, the solitons are a part of the 
general disorder and therefore not visible in any way. The 
change of the structures becomes visible in the angular 
distributions shown in fig. 1111 

Just like in the case without misfit, a further increase of 
the shear velocity leads to a less pronounced effect of the 
shearing again (see curves for v = 1.0 up to u = 3.5) and 



the crystalline structure resumes a greater order. Still, 
due to the shear there is enough disorder in the rows close 
to the walls that there are no distinguishable solitons. 

Here, it is interesting to study the case of very high 
shear velocities (v = 50.0) as the shearing seems to stabi- 
lize the structure with 30 rows. So, if the starting config- 
uration contains 30 rows which contain the same number 
of particles in each row and therefore no solitons, this 
structure is stabilized by the fast moving walls. 

If the starting configuration contains 29 rows and the 
appropriate number of extra particles in the inner rows, a 
different effect can be observed: As the structured walls 
are moving very fast, the rows directly adjacent to the 
walls only feel an averaged-out wall potential and are 
thus not stabilized in a structure with the same num- 
ber of particles as the walls any longer. Therefore it 
becomes energetically preferrable for these rows to con- 
tain the same number of particles like the inner rows, 
which contain extra particles. But as the number of par- 
ticles in the crystal is constant, there are no particles to 
fill up the outer rows. Due to this, holes are created in 
the outer rows which take the role of additional particles. 
These holes can even diffuse further inside the crystal. A 
snapshot showing this effect is presented in fig. [T2] 






FIG. 12. Snapshot of a section of a system with sheared 
walls. The starting configuration contained 29 rows and the 
simulation was carried out at a temperature T — 0.3 and 
a misfit of A = 2.2. The high shear velocity of the walls 
(u = 50.0) leads to the creation of holes. 

Re-running all of these simulations at the same misfit 
and a higher temperature (T ~ 0.5) leads to results that 
can be expected from the case without misfit: As there is 
generally more disorder in the system if the temperature 
is higher, the velocity profiles become more pronounced 
for higher shear rates. In addition, the crystal starts to 
split up into several domains and to change its orientation 
at lower shear rates already than in the case of T = 0.3. 

Another interesting property of these systems is 
the time-dependence of the mean square displacement 
(MSD) of the particles. Here the MSD is defined as 
MSD ~< [xi{t + to) — a;i(io)]^ >, where the average 
< ... > includes an average over all particles in a row 
and over the initial time io- We are only presenting re- 
sults for the MSD in the x-direction parallel to the walls 
here. 

As described above, in cases with a misfit of A = 2.2 
and at low shear velocities {v = 0.1), only the outermost 
row of particles moves along with the walls. This leads to 
an MSD as shown in fig. [T^l where the MSD of the parti- 
cles in the first row is an almost perfect quadratic curve 
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FIG. 13. MSD (in x-direction) of the particles in different 
rows. The starting configuration contained 29 rows in every 
case and the simulation was carried out at a temperature 
T — 0.3 and a misfit of A = 2.2. (a) The shear velocity of 
the walls were chosen as v = 0.1. The oscillation is due to the 
movement of the wall particles which leads to an oscillating 
wall potential for the particles next to the walls, (b) Shear 
velocity v = 0.5, (c) v = 1.2. The inset shows the MSD in the 
case without shearing for comparison. 



of the form MSD = a-t^. As the MSD is the quadratic 
displacement, this means that the fit parameter a — v'^, 
V being the velocity of the particles in this row. As one 
can read off from the graph, the velocity of the particles 
in this case is almost exactly the velocity with which the 
walls are being sheared. The MSD of the rows next to 
the outermost one have significantly smaller values and a 



periodic shape as the particles in these rows are trying to 
cope with the moving walls. The MSD of the inner rows 
is almost constant and independent of the walls, which 
means that the particles just move around randomly. 

At a shear velocity oi v = 0.5 (fig. [T^ middle), where 
the orientation of all the inner parts of the crystal 
changes, the MSD looks completely different: The parti- 
cles in the first row are still moving more than the ones 
in the inner rows, but as the particles of all rows are mix- 
ing in the course of the simulation, the MSDs of different 
rows are getting more similar the longer the simulation 
runs, until finally the particles of all rows are equally 
distributed in the crystal and the MSD of particles stem- 
ming from different rows can not be distinguished any 
more. Besides, all values of the MSDs are very large here 
(note the different scales of the y-axes in fig. [T3)) . 

At even higher shear velocities, where no parts of the 
crystal change their orientations any more, the values of 
the MSD are generally smaller again. As illustrated in 
fig. [T3] (c), the MSD in the row directly adjacent to the 
walls is the second largest one and can not be fitted with a 
quadratic curve of the form MSD = a-t^ any more. This 
is due to the fact that the particles' movement is not pri- 
marily the movement along with the walls any more. The 
MSD in the row next to this row is larger. This is due 
to the incommensurability of the different lattice spac- 
ings. The row directly adjacent to the walls originally 
contained as many particles as the walls while the inner 
rows contain extra particles due to the disappearance of 
one row because of the misfit. In the rows further inside 
of the crystal, the MSD is smaller as the effect of the 
shear and of the incommensurabilities becomes smaller. 

If one compares the MSD of these simulations under 
shear with simulations without shear at the same misfit, 
the effect of the solitons becomes visible (see fig. [131 inset 
of graph c): In the case without shear, the particles in 
the row directly adjacent to the wall do not move much 
as their structure is commensurate with the walls and the 
walls are not moving. Therefore their MSD is the small- 
est one and even shows the overshoot, which one would 
expect to see in a crystal with harmonic potentials where 
the particles are swinging around their ideal lattice po- 
sitions with a certain periodicity. The MSD of the rows 
adjacent to this one are significantly larger in this case as 
the solitons are located in these rows, which causes the 
particles to move around more. In the inner rows, the 
MSD becomes smaller again as there are no solitons in 
these rows. This is clearly different from the values of the 
MSD in the case with shear, where no clearly definable 
solitons are created and therefore all of the outer rows 
exhibit comparable values of the MSD as the incommen- 
surabilities are shared between these rows and lead to a 
lot of movement. Besides this, the shape of the curves is 
different and the absolute values of the MSD are larger as 
the shearing leads to a small, but detectable movement 
along with the sheared walls, while in the case without 
shear no direction is preferred. 

The assignment of particles to rows is done according 



to the initial position of each particle at the beginning of 
the time interval At used as abscissa variable here since 
in the course of time particles change their rows. 



IV. CONCLUDING REMARKS 

In this paper we discussed the influence of shear on the 
structure of a two-dimensional colloidal crystal in con- 
finement between parallel walls. We observed that in the 
case without misfit, at slow shear velocities the layers of 
the crystal are tilted, while the crystal breaks into sev- 
eral pieces with different orientations at medium shear 
velocities and remains almost undisturbed at high shear 
velocities. We showed that the distribution of angles in- 
side of the crystal is an important quantity and gives a 
good insight into the structure of the crystal, along with 
velocity profiles. 

In the case with a misfit, however, the behaviour is 
slightly different, especially at low and at high shear 
rates. At low shear rates, only one single row of particles 
moves along with the walls and at high shear rates holes 
are created in the structure. But the breaking up into 
several pieces also occurs at slighly different shear veloc- 
ities. This is due to the incommensurability in the case 
with a misfit, which would lead to the creation of solitons 



in the case without a misfit. Obviously, the shearing of 
the crystal changes this behaviour significantly. 

Clearly, our results are somewhat qualitative. One 
could in principle run more simulations at different tem- 
peratures and different shear velocities in order to in- 
vestigate the occuring structural transitions (e.g. from 
30 to 29 rows and order-disorder-transitions) with higher 
precision. But this does not seem very sensible at this 
moment because there is no theory with which the re- 
sults could be compared and because experiments with 
such colloidal systems are possible [l-ly, lol-llil l3a |. but 
typically use slightly different parameters for the parti- 
cle interactions, the box sizes and the walls. Therefore 
a rather qualitative study makes sense in this particular 
case and helps to understand the complex interplay of 
confinement and shear and its influence on the structure 
formation of colloidal crystals. 



V. ACKNOWLEDGEMENTS 

This work was performed in the framework of the SFB 
TR6 / project C4 of the Deutsche Forschungsgemein- 
schaft (DFG). It was partially supported by the MAINZ 
graduate school. We are grateful to I. Snook for helpful 
and stimulating discussions. 



[1] K. Zahn, R. Lenke and G. Maret Phys. Rev. Lett. 82 [18 

2721 (1999) 

[2] K. Zahn and G. Maret Phys. Rev. Lett. 85 3656 (2000) [19 

[3] C. Eisenmann, P. Keim, U. Gasser and G. Maret J. 

Phys.: Condens. Matter 16 S4095 (2004) 
[4] K. Zahn, A. Wille, G. Maret, S. Sengupta and P. Nielaba [20 

Phys. Rev. Lett. 90 155506 (2003) 
[5] K. Zahn, J.M. Mendez-AIcaraz and G. Maret Phys. Rev. [21 

Lett. 79 175 (1997) 
[6] H. Konig, R. Hund, K. Zahn and G. Maret Eur. Phys. J. [22 

E 18 287 (2005) 

[7] H. Lowen J. Phys.: Condens. Matter 13 R415 (2001) [23 

[8] C. Likos Phys. Rep. 348 267 (2001) 

[9] M. Koppl, P. Henseler, A. Erbe, P. Nielaba and P. Lei- [24 

derer Phys. Rev. Lett. 97 208302 (2006) 

[10] P. Henseler, A. Erbe, M. Koppl, P. Leiderer and P. [25 

Nielaba Phys. Rev. E 81 041402 (2010) 
[11] L. Assoud, F. Ebert, P. Keim, R. Messina, G. Maret and 

H. Lowen J. Phys.: Condens. Matter 21 464114 (2009) [26 

[12] F. Ebert, P. Dillmann, G. Maret and P. Keim Rev. Set. 

Instr. 80 083902 (2009) [27 

[13] F. Ebert, G. Maret and P. Keim Eur. Phys. J. E 29 311 

(2009) [28 

[14] A. Chowdhury, B.J. Ackerson and N.A. Clark Phys. Rev. 

Lett. 55 833 (1985) [29 

[15] Q.H. Wei, C. Bechinger, D. Rudhardt and P. Leiderer 

Phys. Rev. Lett. 81 2606 (1998) [30 

[16] M. Stieger, P. Lindner and W. Richtering J. Phys.: Con- 
dens. Matter 16 S3861 (2004) [31 
[17] R. Blaak, S. Auer, D. Frenkel and H. Lowen J. Phys.: 
Condens. Matter 16 S3873 (2004) 



R. Blaak, S. Auer, D. Frenkel and H. Lowen Phys. Rev. 
Lett. 93 068303 (2004) 

A. Stipp, R. Biehl, T. Preis, J.N. Liu, A.B. Fontecha, 
H.J. Schope and T. Palberg J. Phys.: Condens. Matter 
16 S3885 (2004) 

M. Stieger, P. Lindner and W. Richtering e-Polymers 046 
(2004) 

G. Petekidis, D. Vlassopoulos and P.N. Pusey J. Phys.: 
Condens. Matter 16 S3955 (2004) 

H. Lowen and C.N. Likos (eds.) J. Phys.: Condens. Mat- 
ter 16 Special Issue (2004) 

A. Ricci, P. Nielaba, S. Sengupta and K. Binder Phys. 
Rev. E 74 010404(R) (2006) 

A. Ricci, P. Nielaba, S. Sengupta and K. Binder Phys. 
Rev. E 75 011405 (2007) 

K. Franzrahe, P. Nielaba, A. Ricci, K. Binder, S. Sen- 
gupta, P. Keim and G. Maret J. Phys. : Condens. Matter 
20 404218 (2008) 

Y-H. Chui, S. Sengupta and K. Binder EPL 83 58004 
(2008) 

Y-H. Chui, S. Sengupta, I.K. Snook and K. Binder Phys. 
Rev. E 81 020403(R) (2010) 

Y-H. Chui, S. Sengupta, I.K. Snook and K. Binder J. 
Chem. Phys. 132 074701 (2010) 

G. Piacente, I.V. Schweigert, J.J. Betouras and F.M. 
Peeters Phys. Rev. B 69 045324 (2004) 
D. Chaudhuri and S. Sengupta Phys. Rev. Lett. 93 
115702 (2004) 

R. Haghgooie and P.S. Doyle Phys. Rev. E 70 061408 
(2004) 



10 



[32] D. Chaudhuri and S. Sengupta J. Chem. Phys. 128 [44 

194702 (2008) 
[33] Y. Han, Y. Shokef, A.M. Alsayed, P. Yunker, T.C. [45 

Lubensky and A.G. Yodh Nature 456 898 (2008) 
[34] C.N. Likos and C.L. Henley Phil Mag. B 68 85 (1993) [46 

[35] L. Assoud, R. Messina and H. Lowen EPL 80 48001 

(2007) 
[36] F. Ebert, P. Keim and G. Maret Eur. Phys. J. E 26 161 [47 

(2008) 
[37] P. Nielaba, K. Binder, D. Chaudhuri, K. Franzrahe, [ 

P. Henseler, M. Lohrer, A. Ricci, S. Sengupta and W. 

Strepp J. Phys.: Condens. Matter 16 S4115 (2004) [49 

[38] K. Franzrahe and P. Nielaba Phys. Rev. E 76 061503 

(2007) [50 

[39] S. Medina Diplomarbeit Johannes-Gutenberg-Universitdt 

Mainz (2010) [51 

[40] S. Medina, P. Virnau and K. Binder J. Phys.: Condens. 

Matter 23 035105 (2011) [52 

[41] A. Imperio and L. Reatto J. Phys.: Condens. Matter 16 

S3769 (2004) 
[42] K. Bagchi, H.C. Andersen and W. Swope Phys. Rev. E [53 

53 3794 (1996) 
[43] HOOMD-blue web page: 

http://codeblue. umich. edu/hoomd-blue\ 



J. A. Anderson, CD. Lorenz and A. Travesset J. Comp. 

Phys. 227 5342 (2008) 

M.P. Allen and D.J. Tildesley Computer Simulation of 

Liquids, Oxford University Press, reprint (1989) 

D.J. Evans and G.P. Morriss Statistical Mechanics of 

Nonequilibrium Liquids, Cambridge University Press, 

2nd edition (2008) 

R. Khare, J. dePablo and A. Yethiraj J. Chem. Phys. 

107 2589 (1997) 

C. Pastorino, K. Binder and M. Miiller Macromolecules 

42 401 (2009) 

C. Pastorino, T. Kreer, M. Miiller and K. Binder Phys. 

Rev. E 76 026706 (2007) 

P.D. Olmsted Rheol Acta 47 283 (2008) and references 

therein 

P. Sollich, F. Lequeux, P. Hebraud and M.E. Gates Phys. 

Rev. Lett. 78 2020 (1997) 

P. Coussot, J.S. Raynaud, F. Bertrand, P. Moucheront, 

J. P. Guilbaud, H.T. Huynh, S. Jarny and D. Lesueur 

Phys. Rev. Lett. 88 218301 (2002) 

A. Sengupta, S. Sengupta and G.I. Menon Phys. Rev. B 

75 180201(R) (2007) and A. Sengupta, S. Sengupta and 

G.I. Menon Phys. Rev. B 81 144521 (2010) 



